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Number theoretic example of scale-free topology inducing self-organized criticality 
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In this work we present a general mechanism by which simple dynamics running on networks 
become self-organized critical for scale free topologies. We illustrate this mechanism with a simple 
arithmetic model of division between integers, the division model. This is the simplest self-organized 
critical model advanced so far, and in this sense it may help to elucidate the mechanism of self- 
organization to criticahty. Its simplicity allows analytical tractability, characterizing several scaling 
relations. Furthermore, its mathematical nature brings about interesting connections between sta- 
tistical physics and number theoretical concepts. We show how this model can be understood as a 
self-organized stochastic process embedded on a network, where the onset of criticality is induced 
by the topology. 

PACS numbers: 89.75.-k, 89. 75. He, 05.65.-l-b, 02.10.De 



In the late 80s Bak, Tang and Wiesenfeld (BTW) [H, Q 
introduced the concept of Self-Organized Criticality 
(SOC) as a mechanism explaining how multicomponent 
systems can evolve naturally into barely stable self- 
organized critical structures without external "tuning" 
of parameters. This single contribution sparkled an 
enormous theoretical and experimental research interest 
in many areas of physics and interdisciplinary science, 
and many natural phenomena were claimed to exhibit 
SOC 3, 4, 5]. However, there was not a general accepted 
definition of what SOC exactly is, and the conditions 
under which it is expected to arise. In order to disengage 
the mechanism of self-organization to criticality one 
should likely focus on rather 'simple' models, and in 
this sense Flyvbjerg recently introduced the "simplest 
SOC model" along with a workable definition of the 
phenomenon @, 0|j namely 'a driven, dissipative system 
consisting of a medium through which disturbances can 
propagate causing a modification of the medium, such 
that eventually, the disturbances are critical, and the 
medium is modified no more —in the statistical sense'. 
On the other hand, in the last years it has been realized 
that the dynamics of processes taking place on networks 
evidence a strong dependence on the network's topology 
[1, Q . Concretely, there exist a current interest on the 
possible relations between SOC behavior and scale-free 
networks Q, characterized by power law degree dis- 
tributions P{k) ~ k~'^, and how self-organized critical 
states can emerge when coupling topology and dynamics 

In this work we introduce a rather simple and general 
mechanism by which the onset of criticality in the 
dynamics of self-organized systems is induced by the 
scale-free topology of the underlying network of interac- 
tions. To illustrate this mechanism we present a simple 
model, the division model from now on, based uniquely 
in the division between integers. We show that this 
model compliances with Flyvbjerg's definition of SOC 



and to our knowledge, constitutes the simplest SOC 
model advanced so far that is also analytically solvable. 
Interestingly, this model establishes connections between 
statistical physics and number theory (see jl4l] for a 
complete bibliography on this topic). 

In number theory, a primitive set of N integers is the 
one for which none of the set elements divide exactly any 
other element jH, [H, Consider an ordered set of 

M — 1 integers {2,3,4, ..,M} (notice that zero and one 
are excluded, and that integers are not repeated), that 
we will name as the pool from now on. Suppose that 
we have extracted N elements from the pool to form 
a primitive set. The division model proceeds then by 
drawing integers at random from the remaining elements 
of the pool and introducing them in the set. Suppose 
that at time t the primitive set contains N{t) elements. 
The algorithm updating rules are the following: 
(Rl) Perturbation: an integer a is drawn from the pool 
at random and introduced in the primitive set. 
(R2) Dissipation: if a divides and/or is divided by say 
s elements of the primitive set, then we say that an 
instantaneous division-avalanche of size s takes place, 
and these latter elements are returned to the pool, 
such that the set remains primitive but with a new size 
N{t+1) = N{t) + l~s. 

This process is then iterated, and we expect the prim- 
itive set to vary in size and composition accordingly. 
The system is driven and dissipative since integers are 
constantly introduced and removed from it, its size 
temporal evolution being characterized by N{t). 

In order to unveil the dynamics undergoing in the model, 
we have performed several Monte Carlo simulations for 
different values of the pool size M. In the upper part 
of fig[T] we have represented for illustration purposes a 
concrete realization of N{t) for M = 10* and A''(0) = 0. 
Note that after a transient, N{t) self-organizes around 
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an average stable value Nc-, fluctuating around it. In the 
inner part of the bottom fig. [H we have plotted in log-log 
the power spectrum of N{t): the system evidences 
noise, with (3 = 1.80 ± 0.01. The former fluctuations 
are indeed related to the fact that at each time step a 
new integer extracted from the pool enters the primitive 
set (external driving Rl). Eventually (according to rule 
R2), a division-avalanche can propagate and cause a 
modification in the size and composition of the primitive 
set. These avalanches constitute the disturbances of the 
system. In figH] (up) we have represented an example 
of the avalanche's size evolution in time. In the same 
figure (bottom) we show the probability P{s) that a 
division- avalanche of size s takes place, for different 
pool sizes M . These latter distributions are power laws 
P{s) ~ s""^ exp(s/so) with r = 2.0 ± 0.1: disturbances 
are thus critical. Observe that the power law relation 
suffers from a crossover to exponential decay at a 
cut-off value Sg due to finite size effects (pool is finite), 
and that the location of these cut-offs scales with the 
system's characteristic size Sq ~ (M/logM)" with 
u) = 1.066 ± 0.003, what is typically characteristic of 
a finite size critical state Q (this characteristic size 
will be explained later in the text). We can conclude 
that according to Flyvbjerg's definition |6], the division 
model exhibits SOC. Division-avalanches lead the system 
to different marginally stable states, that are nothing 
but primitive sets of different sizes and composition. 
Accordingly, for a given pool [2,Af], these time fluctu- 
ations generate a stochastic search in the configuration 
space of primitive sets. 
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FIG. 1: Upper figure: Single realization of the division model 
showing the time evolution of the primitive set size N{t) for 
a pool size M = 10* and A''(0) — 0. Notice that after a 
transient, N{t) self-organizes around an average stable value 
Nc, fluctuating around it. Bottom: (black dots) Scaling be- 
havior of the average stable value Nc as a function of the 
system's characteristic size M/logM. The best fitting pro- 
vides Nc ~ (M/logM)T, with 7 = 1.05 ±0.01. (squares) 
Scaling of Nc as predicted by equation [S] Inner figure: plot 
in log-log of the power spectrum of N{t), showing noise 
with P = 1.80 ± 0.01 (this latter value is the average of 10^ 
realizations of N(t) for 4096 time steps after the transient and 
M = 10''). 



In what follows we discuss analytical insights of the prob- 
lem. Consider the divisor function [20| that provides the 
number of divisors of n, excluding integers 1 and n: 



d{n) = 



k=2 



n-1 



(1) 



where [ J stands for the integer part function. The av- 
erage number of divisors of a given integer in the pool 
[2,M] is then: 
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Accordingly, the mean probability that two numbers a 
and b taken at random from [2,Af] are divisible is ap- 
proximately P = Pr{a\b) + Pr{b\a) ~ 21ogM/M. More- 
over, if we assume that the N elements of the primi- 
tive set are uncorrelated, the probability that a new in- 
teger generates a division-avalanche of size s is on average 
(21ogAf/M)A^. We can consequently build a mean field 
equation for the system's evolution, describing that at 
each time step an integer is introduced in the primitive 



set and a division-avalanche of mean size (2 log M/M)7V 
takes place: 



N{t+l)^N{t) + l 



21ogM 
M 



Nit), (3) 



whose fixed point Nc ~ A// (2 log ill), the stable value 
around which the system self-organizes, scales with the 
system's size as 



Nc{M) 



M 
log M ' 



(4) 



Hitherto, we can conclude that the system's character- 
istic size is not M (pool size) as one should expect in 
the first place, but M/\ogM. This scaling behavior has 
already been noticed in other number-theoretic models 
evidencing collective phenomena [l^, [l^. In figH] we 
have plotted (black dots) the values of Nc as a function 
of the characteristic size M / log M provided by Monte 
Carlo simulations of the model for different pool sizes 
M — 2*, 2^, 2^^ (Nc has been estimated averaging 
N{t) in the steady state). Note that the scaling relation 
|4] holds, however the exact numerical values Nc{M) are 
underestimated by eqlSl This is reasonable since we have 
assumed that the primitive set elements are uncorrelated, 
what is obviously not the case: observe for instance that 
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any prime number p ^ L-^/2j introduced in the primitive 
set will remain there forever. Fortunately this drawback 
of our mean field approximation can be improved by con- 
sidering the function D{n) that defines the exact number 
of divisors that a given integer n e [2, M] has, i.e. the 
amount of numbers in the pool that divide or are divided 
by n: 



D{n) = d{n) + 



M 

n 



- 1. 



(5) 



Define p„ (t) as the probability that the integer n belongs 
at time t to the primitive set. Then, we have 

(6) 

that leads to a stationary survival probability in the 
primitive set: 



1 



f + D{n) 



(7) 



In Fig[3] (right) we depict the stationary survival proba- 
bility of integer n (black dots) obtained through numeri- 
cal simulations for a system with M = 50, while squares 
represent the values of p* as obtained from the eqlTl Note 
that there exists a remarkable agreement. We now can 
proceed to estimate the critical size values Nc{M) as: 



N,{M) 



n=2 



M 



n=2 



1 

D{n) ' 



(8) 



s(t) 



1000 
750 
500 
250 



10000 



20000 t 30000 



40000 



50000 




FIG. 2: Upper figure: Single realization of the division model 
showing the time distribution of division-avalanches. Bottom 
figure: Probability distribution P{s) that a division-avalanche 
of size s takes place in the system, for different pool sizes 
M = 2^° (triangles), M = 2" (inverted triangles), M = 2^^ 
(diamonds) and M = 2^^ (circles). In every case we find 
P(s) ~ s~'^ exp(s/so) with r — 2.0±0.1. Note that the power 
law relation evidences an exponential cut-off due to finite size 
effects at particular values of sq. Inner figure: Scaling of the 
cut-off value so as a function of the system's characteristic 
size M/log Af, with an exponent uj — 1.066 ± 0.003. 



In figH] we have represented (squares) the values of 
Nc{M) predicted by eqlH showing good agreement with 
the numerics (black dots). 

Finally, previous calculations point out that system's 
fluctuations, i.e. division-avalanches distribution P(s) is 
proportional to the percentage of integers having s divi- 
sors. In order to test this conjecture, in figl3] (left) we 
have plotted a histogram describing the amount of in- 
tegers having a given number of divisors, obtained from 
computation of D{n) for M ~ 10®. The tail of this his- 
togram follows a power law with exponent t — 2.0. This 
can be proved analytically as it follows: the responsible 
for the tail of the preceding histogram are those numbers 
that divide many others, i.e. rather small ones (n ^ M). 
A small number n divides typically D{n) ~ [—\- Now, 
how many 'small numbers' have D(n) divisors? The an- 
swer is n, n + 1,..., n + z where 
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The maximum value of z fulfills — M. = i 

n—z n ' 

that is z ~ n^/M. The frequency of D{n) is thus 
fr{D{n)) = n'^/M, but since s = D(n) ~ M/n, we get 



fr{s) ~ Ms ^, and finally normalizing, P{s) ~ s ^. 

Coming back to the Flyvbjerg's definition of SOC, which 
is the medium in the division model? Observe that the 
process can be understood as embedded in a network, 
where nodes are integers, and two nodes are linked if 
they are exactly divisible. The primitive set hence consti- 
tutes a subset of this network, that is dynamically mod- 
ified according to the algorithm's rules. The degree of 
node n is D{n), and consequently the degree distribu- 
tion P{k) ^ k^"^ is scale-free. Hence the SOC behavior, 
which arises due to the divisibility properties of integers, 
can be understood as a sort of anti-percolation process 
taking place in this scale-free network. Observe that the 
division model is a particular case of a more general class 
of self-organized models: a network with M nodes having 
two possible states {on/ off) where the following dynam- 
ics runs: (Rl) perturbation: at each time step a node in 
the state off is randomly chosen and switched on, (R2) 
dissipation: the s neighbors of the perturbed node that 
were in the state on in that time step are switched off, 
and we say that an instantaneous avalanche of size s has 
taken place. Nit) measures the number of nodes in the 
state on as a function of time. Its evolution follows a 
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FIG. 3: Left: Histogram of the amount of integers in [2, 10^] 
that have D divisors. The histogram have been smoothed 
(binned) to reduce scatter. The best fitting provides a power 
law P{D) ~ D"^ with r = 2.01 ± 0.01, in agreement with 
P(s) (see the text). Right: (black dots) Stationary survival 
probability of integer n in a primitive set for a pool size AI = 
50, obtained from Monte Carlo simulations of the model over 
10*' time steps (a preliminary transient of 10'* time steps was 
discarded), (squares) Theoretical prediction of these survival 
probabilities according to equation [T] 



mean field equation that generalizes eq. [3| 

Nit + l) = Nit) + l-^-^N{t), (10) 

where (fc) is the network's mean degree. Accordingly, 
in every case N{t) will self-organize around an average 
value Nc{M). Within regular or random networks, fluc- 
tuations (avalanches) around Nc{M) will follow a Bi- 



nomial or Poisson distribution respectively. However, 
when the network is scale free with degree distribution 
P{k) ^ k ' , fluctuations will follow a power law distri- 
bution P(s) ~ with r = 7, and the dynamics will 
consequently be SOC. In this sense, we claim that scale- 
free topology induces criticality. 

Some questions concerning this new mechanism can be 
depicted, namely: which is the relation between the spe- 
cific topology of scale-free networks and the power spec- 
tra of the system's dynamics? Which physical or natural 
systems evidence this behavior? 

With regard to the division model, the bridge between 
statistical physics and number theory should also be 
investigated in depth. This includes possible general- 
izations of this model to other related sets such as k- 
primitive sets [2l| . where every number divides or is di- 
vided by at least k others {k acting as a threshold pa- 
rameter), to relatively primitive sets ^25| and to cross- 
primitive sets (where this will introduce coupled SOC 
models). From the computational viewpoint [23| . proper- 
ties of the model as a primitive set generator should also 
be studied. Of special interest is the task of determining 
the maximal size of a fc-primitive set [l6l . [2l| , something 
that can be studied within the division model through 
extreme value theory [4]. 
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